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ABSTRACT 

We present a very large high-resolution cosmological N-body simulation, the 
Millennium-XXL or MXXL, which uses 303 billion particles to represent the forma- 
tion of dark matter structures throughout a 4.1 Gpc box in a ACDM cosmology. We 
create sky maps and identify large samples of galaxy clusters using surrogates for 
four different observables: richness estimated from galaxy surveys, X-ray luminosity, 
integrated Sunyaev-Zeldovich signal, and lensing mass. The unprecedented combina- 
tion of volume and resolution allows us to explore in detail how these observables 
scale with each other and with cluster mass. The scatter correlates between different 
mass-observable relations because of common sensitivities to the internal structure, 
orientation and environment of clusters, as well as to line-of-sight superposition of 
uncorrelatcd structure. We show that this can account for the apparent discrepancies 
uncovered recently between the mean thermal SZ signals measured for optically and 
X-ray selected clusters by stacking data from the Planck satellite. Related systematics 
can also affect inferences from extreme clusters detected at high redshift. Our results 
illustrate that cosmological conclusions from galaxy cluster surveys depend critically 
on proper modelling, not only of the relevant physics, but also of the full distribu- 
tion of the observables and of the selection biases induced by cluster identification 
procedures. 
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1 INTRODUCTION 

Confrontation of observational data with theoretical models, 
and in particular with numerical simulations, has been a key 
factor enabling the rapid recent progress of cosmological re- 
search. Without it, we would have not arrived at the current 
structure formation paradigm, which is now being subjected 
to ever more detailed scrutiny. Further development of this 
fruitful approach requires current and future observations 
to be matched with equally precise theoretical models. The 
order of magnitude advances made by new surveys hence 
require new simulations with comparable improvements in 
statistical power and accuracy. 

An inevitable consequence of the increasing accuracy of 
observational data and the growing sophistication of numeri- 
cal simulations is that comparing them becomes a non-trivial 
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task in its own right. It has long been appreciated that the 
distribution of properties in a sample of observed objects is 
shaped not only by the relevant physics but also by the ob- 
servational methods used to detect and characterise them. 
The resulting measurement biases have often been neglected 
in the past, but this is no longer possible in the era of 'preci- 
sion cosmology' where the systematic errors in observational 
results are typically comparable to or larger than their sta- 
tistical errors. Detailed modelling of a given observational 
programme is not optional in this situation, but rather is an 
unavoidable step in the proper interpretation and exploita- 
tion of the data. 

In this paper, we present a major new effort in this 
direction, aiming to address two aspects of the physics 
of galaxy clusters that have recently attracted a lot of 
interest. The first concerns the relations between opti- 
cal richness, lensing mass, X-ray luminosity and thermal 
Sunyaev-Zeldovich (tSZ) decrement. It is critically impor- 
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tant to calibrate how these observables scale with 'true' 
mass if cluster counts are to be used to place robust con- 
straints on cosmological parameters. The Planck Collab- 
oration has recently reported puzzling inconsistencies in 
the scaling relations measured for different samples, sug- 
gesting an unexp ected dichotomy in the gas prop erties 
of galaxy clusters (jPlanck Collaboration et al.ll2011q ). The 
other aspect concerns the inferred masses of extreme galaxy 
clusters. This is interesting because discoveries of mas- 
sive clusters at high redshift have repeatedly been sug- 
gested to be in tension with the standard ACDM model, 
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To study these questions, we use a new state-of-the-art 
simulation of the evolution of the dark matter structure that 
provides the arena for the formation and evolution of galax- 
ies. This Millennium-XXL is the largest high-resolution cos- 
mological N-body simulation to date, extending and comple- 
menti ng the previous Millennium and Millennium -II simula- 
tions (|Springelll20"05l ; iBovlan-Kolchin et aLlfeood) . It follows 
the dark matter distribution throughout a volume equiv- 
alent to that of the whole sky up to redshift z = 0.7, or 
equivalently, of an octant up to redshift z — 1.4. Its time 
and mass resolution are high enough to allow detailed mod- 
elling of the formation of the galaxy populations targeted 
by future large surveys, as well as of the internal structure 
of extremely rare and massive clusters. It is also well suited 
to studying a number of other probes of the cosmic expan- 
sion and structural growth histories, for example, baryonic 
acoustic oscillations (BAOs), redshift space distortions, clus- 
ter number counts, weak gravitational lensing, and the inte- 
grated Sachs- Wolfe effect. The volume of Millennium-XXL 
is very much larger than can be followed by direct hydro- 
dynamical computations, but its resolution is sufficient for 
galaxy formation to be followed in detail within each halo 
by applying semi-analytic models to its merger tree. 

Using mock observations of galaxy clusters in the 
Millennium-XXL, we show in this study that current in- 
terpretations of cluster surveys are significantly affected by 
systematic biases. In particular, we show that the apparent 
inconsistencies highlighted by the Planck Collaboration in 
the mean SZ and X-ray signals measured for optically and 
X-ray selected cluster samples can be understood as result- 
ing from substantial and correlated scatter in the various 
observables among clusters of given 'true' mass. Currently 
there appears to be no compelling evidence for unknown 
processes affecting the gas properties of clusters or for a bi- 
modality in cluster scaling properties. We also comment on 
the implications of our results for constraining cosmology 
using extreme clusters at high redshift. 

Our paper is structured as follows. Section 2 is devoted 
to presenting the Millennium-XXL and our techniques for 
modelling the observable properties of galaxy clusters. In 
particular, Section 2.1 provides technical and numerical de- 
tails of the simulation, while Section 2.2 describes our sur- 
rogates for X-ray luminosity, gravitational lensing mass, op- 
tical richness and thermal SZ flux. In Section 3, we then 
explore the impact and implications of various selection bi- 
ases. We explain how we identify clusters in Section 3.1, and 
in Section 3.2 we analyse the bulk of the cluster population, 



with a focus on extreme objects. We discuss the implication 
of our findings for the conundrum reported by the Planck 
Collaboration in Section 4. Finally, we present our conclu- 
sions in Section 5. 



2 NUMERICAL METHODS 

In this section we describe our dark matter simulation and 
the way we use it to construct surrogates for four observa- 
tional properties of galaxy clusters; their optical richness, 
their weak gravitational lensing signal, their X-ray luminos- 
ity and their thermal Sunyaev-Zeldovich (tSZ) amplitude. 



2.1 The MXXL N-body simulation 

The 'Millennium-XXL Simulation' (MXXL) follows the non- 
linear growth of dark matter structure within a cubic region 
of 4.11 Gpc (3/i _1 Gpc) on a side. The dark matter distri- 
bution is represented by 6720 3 = 303, 464, 448, 000 particles, 
substantially exceedin g the number used in all previous sim- 
ulations of this type fSpringel et all 120051 : iKim et al.ll2009l : 
iTevssier et all 120091 : IPrada et al. 201 1 | ) ap art from the re- 
cent 'Horizon Run 3' of IKim et all (|201lh . which has 20% 
more particles but 40 times poorer mass resolution. We 
note that the simulated volume of the MXXL is equiva- 
lent to that of the whole observable Universe up to red- 
shift z = 0.72. It is more than 200 times that of our 'Mil- 
lennium Simulation' (MS, ISpringel et al.ll2005h . almost 30 
times that of the r ecently completed MultiDark simulation 
l|Prada et al.ll201lf ) but still only 2% that of the Horizon 
Run 3. The MXXL is also about 7 times larger than the ex- 
pected volum e of the Baryon Osc illation Spectroscopic Sur- 
vey (BOSS) (|Schlegel et all 120071 ) and about twice that of 
the planed JPAS^. Its particle mass is m p = 8.456 xlO 9 M , 
approximately 7 times that of the MS but more than 300 
times smaller than that of the 'Hubble Volume Simulation' 
(jEvrard et al.ll2002f ). completed a decade ago with a com- 
parable volume to MXXL. The mass resolution of MXXL 
is sufficient to identify the dark matter haloes hosting cen- 
tral galaxies with ste llar mass exceeding ~ 1.5 x 10 10 Mq 
(|De Lucia et al.ll2006h . and also to predict robustly the in- 
ternal properties of the haloes corresponding to very massive 
clusters, which are represented by more than 100, 000 dark 
matter particles. The Plummer-equivalent softening length 
of the gravitational force is e = 13.7 kpc, which translates 
into a dynamic range of 300, 000 per dimension, or formally 
to more than 2 x 10 16 resolution elements within the full 
simulation volume. This large dynamic range can be appre- 
ciated in Fig. [1] where we show the large-scale density field 
together with the internal structure of a few selected massive 
clusters. 

The MXXL adopts a ACDM cosmology with the 
same cosmological parameters and output times as the 
previous two Millenn i um si mulations l|Springel et al.ll2005l : 
IBovlan-Kolchin et al.ll200d l. This facilitates the joint use 
of all these simulations in building models for the galaxy 
population. Specifically, the total matter density, in units 
of the critical density, is fi m = Qdm + Sib = 0.25, where 
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Qb = 0.045 refers to baryons (although these are not ex- 
plicitly treated in the simulation); a cosmological constant, 
SIa = 0.75, gives a flat space geometry; the rms linear den- 
sity fluctuation in 10.96 Mpc spheres, extrapolated to the 
present epoch, is erg = 0.9; and the present-day Hubble con- 
stant is Ho — 73 km s~ Mpc - . Although this set of param- 
eters is discrepant at about the 3cr level wit h the latest con- 
strain ts from CMB and LSS observation s (iKomatsu et al. 
l201lT) . the scal ing technique pr oposed by lAngulo fc White 
(j2O10T) (see also lRuiz et alj|201ll ) allows the Millennium sim- 
ulations to provide theoretical models for the formation, evo- 
lution and clustering of galaxies over the full range of cos- 
mologies allowed by current observational constraints. The 
parameter offset with respect to the best current observa- 
tional estimates lies mainly in the high value for erg, but 
this is an advantage for reliable scaling of the simulation re- 
sults to other cosmologies, as this requires interpolation on 
the stored MXXL/MS/MS-II data which is only possible for 
target cosmologies with lower <rg than used in the MS. 

2.1.1 Initial conditions 

The initial unperturbed particle load for the simulation was 
built by periodically replicating a 280 3 particle cubic glass 
file twenty-four times in each coordinate direction. T he glass 
file wa s created for the MXXL usin g the method of I White] 
l |l996h (see also iBaugh et al.lll995h . The initial conditions 
were then produced by computing displacements and veloc- 
ities for each of the particles at starting redshift z s tart = 63, 
using an upgraded versio n of the code origina lly developed 
for the Aquarius Project (|Springel et al.ll2008l ). Further im- 
provements include communication and memory optimisa- 
tions, as well as the us e of second-order L agrangian pertur- 
bation theory (2LPT, [g coccimarr ol ll998h . rather than the 
Zel'dovich approximation for computing the position and 
velocity perturbations of each particle. The latter modifi- 
cation is particularly important for the present study, since 
the abundance of high mass haloes is sensitive to initial tran- 
sients, which a re much smaller and decay more quickly when 
2LPT is used (|Crocce et al.ll2006h . 

Another important change is the introduction of a new 
approach to generate Gaussian initial fluctuations. Rather 
than setting the phases of the modes in fc-space (as done, 
for example, in the MS), we first generated a real-space 
white noise field. The Fourier transform of this field was 
then used to set the amplit udes of all of the modes needed to 
make the initial conditions (|Salmonlll996l ; iBertschingeihoOll ; 
lHahn fc Abel 1201 if ). For the MXXL, the real-space white 
noise field was created on a 9216 3 grid. Only modes within 
a spherical fe-space volume of radius 6720/2 = 3360 times 
the fundamental frequency (i.e. below the particle Nyquist 
frequency) were used to generate the displacement and ve- 
locity fields (all other modes were given zero amplitude). 

The use of a white noise field in real space, while not 
necessary for the MXXL initial conditions themselves, will 
make it much easier to resimulate arbitrary MXXL regions 
of interest at higher resolution, for example, the extreme 
objects illustrated in the present paper. This is because in 
our new approach it is unnecessary to reproduce the en- 
tire white noise field at the original resolution in order to 
capture the phases of large-scale modes. A consequence is 
the ability to create consistent sets of initial conditions for 



resimulations (including 'resimulations of resimulations' at 
yet higher resolution) for arbitrarily defined subregions over 
a huge dynamic range. The real-space white noise field is 
generated in a special top-down hierarchical fashion, based 
on an oct-tree, making it easy to generate coarse represen- 
tations of the MXXL field at low computational cost. The 
MXXL white noise field itself occupies just a small subvol- 
ume of a single realisation of a huge white noise field created 
in a hierarchical way. This realisation is specified everywhere 
to a resolution below the likely free streaming scale of cold 
dark matter. This means that resimulations of parts of the 
MXXL volume can be created at any desired resolution as 
the phases are fully specified everywhere in advance. A full 
description of this method will be given in Jenkins (2012, in 
preparation). 

2.1.2 The simulation code 

Evolving the distribution of the dark matter particles in 
the MXXL under their mutual gravitational influence was 
a formidable computational problem. Storing the positions 
and velocities of the particles in single precision already re- 
quires about 7 TB of memory. As each particle exerts a force 
on every other particle, a CPU- and memory-efficient ap- 
proximate calculation of the forces is of paramount impor- 
tance. It is also necessary to develop new strategies to deal 
with the huge data volume produced by the simulation. Us- 
ing the same analysis approach as for the Millennium Sim- 
ulation would have resulted in more than 700 TB of data, 
adding a severe data analysis problem and significant disk 
space costs to the computational challenge. 

In order to alleviate these problems, we developed a spe- 
cial "lean" version of the Tree-PM code GADGET-3, which 
improves the scalability and memory efficiency of the code 
considerably, outperforming the highly optimised version of 
GADGET-2 (|Springell [2003 ) used for the MS. GADGET-3 
computes gravitational forces with a TreePM method by 
combining a particle-mesh (PM) scheme with a hierarchi- 
cal tree-method, and it uses spatially and temporally adap- 
tive time-stepping, so that short time-steps are used only 
when particles enter localised dense regions where dynam- 
ical times are short. A significant improvement in the new 
code is a domain decomposition that produces almost ideal 
scaling on massively parallel computers. Finally, the MXXL 
version of GADGET-3 uses aggressive strategies to minimise 
memory consumption without compromising integration ac- 
curacy and computational speed. 

The code also carries out a significant part of the re- 
quired post-processing on-the-ny as an integral part of the 
simulati on. This includes group-finding via the Friends-of- 
Friends (|Davis et al.lll985T) (FOF) algorithm, application of 
the SUBFIND algorithm |Springel et al.ll200lT ) to find gravi- 
tationally bound subhaloes within these groups, and calcu- 
lation of basic properties of these (sub)haloes, like maximum 
circular velocities, cumulative density profiles, halo shapes 
and orientations, velocity dispersions, etc. These extended 
halo and subhalo catalogues are then stored at the same out- 
put times as for the other Millennium simulations, allowing 
the construction of detailed (sub)halo merger trees. Full par- 
ticle data are, however, stored only at a handful of redshifts, 
very significantly reducing the stored data volume. Only 72 
bytes per particle are needed by the simulation code during 
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Figure 1. The projected density of dark matter in the MXXL simulation at z = 0.25. The insets correspond to circles of radius 5.5 Mpc 
centred on the most extreme clusters identified according to our surrogates for X-ray luminosity, optical richness, lensing signal and 
integrated thermal tSZ strength (see section [2.3l l. The underlying image is a projection of the dark matter density in a slab of thickness 
27 Mpc, and width 2050 Mpc. It is oriented so that it contains three of the selected clusters, as indicated in the figure (the lensing 
example is from a different slice). The whole simulation box is actually twice as wide, spanning 4110 Mpc. All four cluster images and 
the large-scale slice use the same colour scale, which varies in shade from light blue in the least dense regions to orange and white in the 
densest regions. 

Also, this makes it easier for our code to reach close-to- 
optimum work- and load-balance during the calculation. 

The final production run carried out approximately 87 
trillion force calculations to reach 2 = 0, and used about 
28.5 TB of RAM, nearly the whole available physical mem- 
ory of JuRoPa. The run-time was 9.3 days (wall-clock), 
equivalent to 2.86 million CPU hours (or 326 years) in serial. 
Of this time, 15% were required for running our on-the-fly 
postprocessing software, notably the group finding, the sub- 
structure finding, and the power spectrum calculation, and 
another 14% were needed for I/O operations. The total long- 
term storage space required for all MXXL data products is 
about 100 TB, down by a factor of about 8 per particle rela- 
tive to the approach used for the MS and MS-II simulations. 

2.2 Basic validation results 

At redshift z — 0, the MXXL contains more than 700 million 
haloes with at least 20 particles. These account for 44% of all 
the mass in the simulation. Among these objects, 23 million 



normal dynamical evolution. When the in-lined group and 
substructure finders are enabled as well (which is optional) , 
the peak memory consumption per particle increases by a 
further 26 bytes. 



2.1.3 Computational cost and code performance 

The MXXL simulation was carried out in the late summer 
of 2010 on the JuRoPa machine at the Jiilich Supercom- 
puting Centre (JSC) in Germany. A partition of 1536 com- 
pute nodes was used, each equipped with two quad-core Intel 
X5570 processors and 24 GB of RAM. We ran our code in 
a hybrid MPI/shared memory setup on 12888 cores, placing 
one MPI task per processor socket (3072 in total), and em- 
ploying all four cores of each socket via threads. This setup 
turned out to be advantageous compared with a pure MPI- 
parallelisation based on 12888 MPI tasks, because it reduces 
the amount of intra-node MPI communication, and min- 
imises the RAM required for MPI communication buffers. 



Figure 2. The differential Friends-of-Friends (FOF) halo mass 
function (top panel) of the MXXL (blue), MS (red) and MS-II 
(green). The MXXL provides vastly superior sampling of the mas- 
sive end, where the abundance of objects drops exponentially as a 
function of mass. Combined, the three simulations cover about 8 
decades in halo mass. The vertical lines mark the halo resolution 
limits (20 particles) of the three simulations. For comparison, we 
also display a fit to the mass function of all self-bound subhaloes 
in the three Millennium simulations (dashed). The bottom panel 
gives the ratio of the three mass functions to an analytic fitting 
formula given in the text. We see that the simulations agree ac- 
curately with each other for intermediate masses, but also that 
different methods for identifying structures disagree significantly 
in the expected number density of objects of given mass, espe- 
cially at the high-mass end. 



have a value of M20JEI larger than that of the Milky- Way's 
halo (M200 = 2 x 10 12 Mq) and 464 have a value in excess 
of that of the Coma galaxy cluster (M200 = 2 x 10 15 M Q ). 
In Fig. [2] we show the differential halo mass function (FoF 
masses for a linking length b = 0.2) at the present epoch, 
which is a robust way of describing the abundance o f non- 
linear objects as a function of mass (|Davis et al.ll 19851 ). The 
most massive halo at z = has Mf f = 8.98 x 10 15 Mq . Such 
extreme objects are so rare that they can only be found in 
volumes as large as that of the MXXL. We compare the 
MXXL results with similar measurements from the MS and 
MS-II simulations. For masses where the three simulations 
have good statistics and are away from their resolution lim- 
its, the agreement is at the few percent level. The results 
from all three simulations are well described by 
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2 We define the conventional virial mass of a halo M200 as the 
mass within a sphere centred on the potential minimum which 
has mean density 200 times the critical value. 




Figure 3. Projected dark matter density for the 15 most mas- 
sive MXXL haloes (according to A/200 at 2 = 0.25. Each image 
corresponds to a region of dimensions 6 X 3.7 /i — x Mpc wide and 
20 h~ 1 Mpc deep. Note the large variation in shape and internal 
structure among these clusters. In particular, the most massive 
cluster, shown in the top left corner, has no clear centre but rather 
displays several distinct density peaks of similar amplitude. 



where po is the mean mass density of the universe, a(M) is 
the variance of the linear density field within a top-hat filter 
containing mass M, and f(a) is the fitting function 
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The residuals from this analytic halo mass function, dis- 
played in the bottom panel of Fig. [2] show that it de- 
scribes the numerical results accurately (to better than 5% 
over most of the mass range) over eight orders of magni- 
tude in halo mass, extending the accur acy of previous mod- 
els to larger and to smaller scales (e.g . I Jenkins et al.ll200ll ; 
IWarren et aT]|2006l ; PrTnker et al.ll2008l ). In Fig.[2]the dashed 
lines show a fit of this same analytic form to the mass func- 
tion of all self-bound subhaloes (as identified by SUBFIND) 
in the MXXL, MS and MS-II simulations. These curves cor- 
respond to the fit 
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The difference between the two fits illustrates how the mass 
function of objects depends on the way they are defined. 
This is especially important at the high-mass end. For ex- 
ample, the expected abundance of haloes with M ~ 10 15 Mq 
changes by a factor of ~ 2 when FoF haloes and self-bound 
subhaloes are compared. 

The difficulty in unambiguously defining haloes and 
their associated mass is in part a consequence of the fact 
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Figure 4. Matter power spectra measured directly in the MXXL, 
MS, and MS-II top panel). The black line shows the power spec- 
trum used to generate the initial conditions, linearly evolved to 
z = 0. The dashed lines show the Poisson power level of each 
simulation, which becomes significant only at the smallest scales. 
The Poisson power has been subtracted from the measured power 
spectra in this figure. In the bottom panel, we show the ratio of 
the measured power spectra to the actual realisation of the linear 
theory used to generate the initial conditions of each simulation. 
This procedure reduces sampling noise due to the finite number of 
modes at small wavenumber. The arrows mark the gravitational 
resolution limits (2tt /softening length) of the three simulations. 



that large haloes do not form a homogeneous population. 
In fact, they display considerable variety in structure and 
environment. We illustrate this in Fig. [3] by showing the 
fifteen most massive clusters in the MXXL at z = 0.25, 
selected according to M200. Among this group there is con- 
siderable diversity in shape, concentration and the amount 
of substructure, despite all the objects having very similar 
virial mass. This already suggests that careful modelling of 
mass estimators will be needed to compare numerical simu- 
lations with observed massive clusters at high redshift. Small 
changes in the estimated mass of an object can dramatically 
change the predicted probability of its existence within any 
given cosmological model. 

The diversity of massive clusters may also have impor- 
tant consequences for other observational studies. Matched 
filters are often applied to data in order to maximise the 
signal-to-nois e ratio of, for instance, weak lensing or tSZ de- 
tectio ns (e.g. ISchneider|[l99r3 ; iMelin et al.1l200fj| ; iRozo et al.l 
2011). Such filters use a model for the spatial distribution of 
the signal as prior information (e.g. in the form of density 
or pressure profiles) but in many cases (and in particular for 
the most massive objects) the structure of individual clus- 
ters will not conform to these assumptions. For instance the 
top left halo in Fig. [3] which is the most massive cluster 
in the MXXL at z = 0.25, does not have a clear centre. 
In such cases, the signal may be seriously misestimated by 



a matched filter, potentially biasing cosmological inferences 
from the measurement. 

In Fig. |3J we show power spectra of the mass density 
field at the present epoch. The results are a combination 
of two measurements. Large-scale modes were computed us- 
ing a global 9216 3 mesh, whereas the mean amplitude of 
smaller modes was calculated by folding the density field 
64 times alo ng each direction an d projecting it onto a new 
9216 3 mesh (| Jenkins et al.lll99ct ). This method effectively 
reaches the same spatial resolution as a 589, 824 3 mesh. For 
comparison, we also show results for the MS and MS-II sim- 
ulations. Clearly, only the MXXL simulation probes scales 
significantly beyond the turnover in the power spectrum. 
The MXXL is also the only one among the three runs that 
provides good sampling of the baryonic acoustic oscillations 
(BAO). We note that at low redshift these features already 
show clear signs of b eing affected by nonlinear evolution (e.g. 
lAngulo et al.l [20081 ) . making the MXXL particularly valu- 
able for studying systematic effects in large-scale galaxy sur- 
veys aiming at precise measurements of the BAO features. 
Throughout the nonlinear regime, the power spectra of the 
three Millennium simulations show excellent agreement up 
to the scales where the spatial resolution limits of each run 
kick in, manifested as a reduction in power relative to higher 
resolution simulations. 



2.3 Surrogate observables 

Using the DM distribution and halo catalogues described 
in the previous section we have created surrogate observ- 
ables that mimic the four main techniques used observa- 
tionally to discover and characterise large clusters: optical 
galaxy counts, gravitational lensing, X-ray emission, and the 
Sunyaev-Zel'dovich signal imprinted on the microwave back- 
ground radiation. 

Rather than attempting to follow the baryonic physics 
directly in the simulation, we have constructed simple prox- 
ies for these observables, based directly on the dark mat- 
ter distribution. This necessarily schematic approach avoids 
the uncertainties of any specific implementation of baryonic 
processes such as star and black hole formation and the as- 
sociated feedback, while allowing us to take advantage of 
the characteristics of the MXXL, namely its combination of 
very large volume and relatively high mass resolution. Our 
approach can easily be updated as a better understanding 
of the relation between the dark matter structure of galaxy 
clusters and any particular observable is achieved. Our main 
goal in this paper is not to produce accurate a priori pre- 
dictions for the observables, but to look for surrogates that 
correctly rank the expected signal strengths and represent 
the scatter and the correlations between observables in a 
realistic way. We can then study the diversity of clusters 
and quantify the extent to which different methods select 
different cluster populations. 

We focus our analysis on redshift z = 0.25 because 
the most massive halo in t he observable Universe should 
be roughly at that redshift |Holz fc Perlmutter|[2O10l ). This 
redshift also corresponds to the median redshift of galax- 
ies in the pho tometric catalogu e of the Sloan Digital Sky 
Survey (SDSS|York et al.|[2000l ). which provides one of the 
largest samp les of optically detecte d clusters - the MaxBCG 
catalogue of iKoester et al.l (|2007bl ), which has been widely 
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used to compute scaling relations for optically selected clus- 
ters. Our results do not change qualitatively if we pick an- 
other redshift between z = and 1. The largest simu- 
lated haloes should resemble the most massive observable 
objects because the MXXL volume is comparable to or ex- 
ceeds that accessible to real surveys. In 1000 all-sky light- 
cones up to z = 0.6 (built by placing observers randomly 
on 'Milky Way' haloes), we find the most massive halo is 
typically between z — 0.1 and z = 0.3, and has a virial mass 
M 2 oo ~4x 10 15 M , roughly 75% that of the most massive 
MXXL ha lo at z = 0.25, consisten t with previous analytic 
estimates (|Holz fc PerlmutteJl20ld ). 

Finally, we note that we normalise our surrogates to 
match observed scaling relations between observables and 
halo mass (which we discuss in Section 4). This allows a 
direct comparison of population properties to observational 
data, side-stepping issues of possible offsets due to incorrect 
cosmological parameters, to the schematic nature of our sur- 
rogates, or to observational details such as filter shapes. We 
now outline how we construct 2D maps from which we can 
identify clusters and measure our various surrogate observ- 
ables. 



2.3.1 Optical maps 

The first observational approach we consider is the de- 
tection of rich clusters in optical surveys, which relies 
on finding large groups of galaxies in a narrow redshift 
range and at similar projected positions on the sky. In 
order to mimic this, we start by constructing a three- 
dimensional galaxy ca talogue using a halo occupation distri- 
bution (HOD) model dKauffmann et al.lll997l ; 15 enson et al.l 
l2000l ; |Peacock fc Smithll2000l ) to populate each MXXL halo 
with galaxies. We assume that every halo with M200 above 
1.4xl0 12 M Q hosts one central galaxy and a number of satel- 
lites drawn from a Poisson distribution with mean equal to 



the halo mass in units of 5.7 x 10 12 Mq raised to the 0.9 
power. This HOD is similar to th at derived for red galaxies 
in the SDSS |Zehavi et alj|201ll ). but it is tuned to repro- 
duce the obser ved mass-richness relat ion for galaxy clusters 
as measured bv ljohnston et all (|2007T l. The central galaxy is 
placed at the minimum of the gravitational potential of the 
dominant SUBFIND substructure in the halo, whereas the 
satellite galaxies are identified with randomly chosen dark 
matter particles of the FoF group. The latter ensures that 
the effects of halo ellipticity, alignment and substructure, 
which are importan t for reproducing the small-scale corre - 
lations of galaxies (|Zu et al.l [20081 ; Ivan Daalen et all l201lh 
are included in our modelling. The resulting catalogue at 
z — 0.25 contains more than 150 million galaxies (50% of 
which are satellites). 

We note that the basic assumption of this HOD mod- 
elling, namely that the galaxy content of a halo is statis- 
tically determined exclusively by its mass, is not expected 
to hold in detail. Models that follow galaxy formation ex- 
plicitly predict a dependence of the HOD on other proper- 
ties such as the halo formation tim e (e.g. iGao et ah I l2005l ; 
IZhu et al.|[2006l ; ICroton et al.ll2007h . Nevertheless, these ef- 
fects are weak so we do not expect them to affect our con- 
clusions. Slightly different HODs change the average number 
of galaxies in our clusters but make little difference to the 



correlation of optical properties with other aspects of the 
halo. 



2.3.2 Lensmg maps 

The second identification approach we consider is weak grav- 
itational lensing. Although direct mass measurements using 
this effect are only possible for the most massive individual 
clusters, lensing can be used to estimate precise mean masses 
by stacking a large number of clusters selected according to 
speci fic criteria (e.g. Mandclb aum et ah 2006; Sheldo n" et alj 
l2009h . It is easy to see that orientation effects and both cor- 
related and uncorrelated large-scale structure can play an 
important role in defining the lensing signal of an individual 
cluster. As a result, large N-body simulations where the full 
line-of-sight density distribution is properly modelled are 
needed to calculate accurately the distribution of expected 
signals. 

Our modelling of lensing maps uses the distant observer 
approximation. We neglect the evolution of clustering along 
the line-of-sight and assume that all the mass along a line-of- 
sight from z = to z = 1.37 (the side length of the MXXL) 
contributes equally to the convergence field. Under these 
assumptions we create "weak lensing mass maps" by pro- 
jecting the simulated mass density of the z — 0.25 snapshot 
along one axis. MXXL particles are mapped onto a 32, 768 2 
mesh using a nearest grid point (NGP) mass assignment 
scheme, yielding an effective transverse spatial resolution of 
~ 92/i _1 kpc. 

A more realistic approach would vary the weight as- 
signed to mass at different redshifts assuming a specific red- 
shift distribution for the background source galaxies. The 
highest weight would go to material which is "halfway" to 
the sources. Our 4.1 Gpc projection length will clearly tend 
to overestimate projection effects from distant matter. How- 
ever, it turns out that structures far in front or far behind 
the clusters produce only a small fraction of the projection 
effects; most come from their immediate surroundings and 
from the directional dependence of the projection of their 
internal structure. Our simple model can be regarded as 
treating these aspects quite accurately and as giving an over- 
estimate of the (subdominant) effects of distant projections. 



2.3.3 X-ray maps 

Another important route to detecting and characterising 
massive clusters is through X-ray emission from their hot 
intracluster gas. The local X-ray emissivity is proportional 
to the square of the gas density and (approximately) to the 
square root of its temperature. This makes X-rays a partic- 
ularly sensitive tracer of the inner regions of clusters, where 
densities can be thousands of times higher than in the out- 
skirts. Simulations with high spatial and mass resolution are 
needed to probe these inner regions adequately. 

We estimate a density local to each particle by using 
kernel-interpolation over its 32 nearest neighbours, a com- 
mon method in SPH calculations, while we take the local 
temperature to be proportional to the velocity dispersion of 
the subhalo in which the particle is located. Particles out- 
side subhaloes are taken to have zero temperature. With 
these quantities in hand, we compute a 32, 768 2 pixel X-ray 
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map by summing up the density times the square root of the 
temperature for all the particles along a given line of sight. 

This surrogate clearly neglects dynamical effects on 
cluster luminosity during vi olent cluster merge rs. Recent hy- 
drodynamical simulations (|Rasia et al.l 1201 lh suggest that 
luminosity enhancements during such events can be sub- 
stantial, particularly for high Mach number (> 2.5) and 
for equal progenitor masses. On the other hand, observa- 
tional data suggest that disturbed, apparently merging clus- 
ters tend to have lower than average X-ray luminosities for 
their mass, whereas symmetric equilibrium clusters often 
have cold cores and thus higher than average luminosities 
jArnaud et alj|2010h . Even the sign of merger effects thus 
seems unclear. 



2.3-4 Thermal SZ maps 

The final cluster property we co nsider is their ther- 
mal S unyae v-Zel'dovich (tSZ) signal |Sunvaev fc Zeldovichl 
1 19721 . 1 19801 ). This effect causes a characteristic distortion 
of the spectral shape of the cosmic microwave background 
(CMB) as a result of inverse Compton scattering of CMB 
photons off the electrons in the hot intracluster plasma. The 
integrated magnitude of the effect is proportional to the to- 
tal thermal energy content of the hot electrons in the cluster, 
or, equivalently, to the gas mass times the mean gas temper- 
ature. Along any given line-of-sight the effect is proportional 
to the line integral of the gas pressure. 

In our analysis, we again assume the gas to be dis- 
tributed like the dark matter and to be isothermal within 
each quasi-equilibrium subhalo, associating a temperature 
to each simulation particle proportional to the velocity dis- 
persion of its host subhalo. We then create a 32, 768 2 pixel 
tSZ map by projecting the thermal energies of all MXXL 
particles along one of the box axes. 



3 RESULTS 

In this section we use our simulated halo catalogues and 
mock observational maps to examine various systematic ef- 
fects that can have an impact on the scatter and mean am- 
plitude of cluster scaling relations. We cumulatively include 
effects due to sample selection, to spurious cluster identifi- 
cation resulting from projection effects, to misidentification 
of cluster centres, and to contamination by structures along 
the line-of-sight, examining in each case the impact on rela- 
tions between mean tSZ and optical richness, and between 
weak lensing mass and optical richness. 



3.1 Cluster catalogues 

We first specify how we identify optical clusters in our 
galaxy map. We have implemented a group finder similar 
to that employed t o build the MaxBCG cluster catalogue 
from SDSS galaxies (|Koester et al.|[2~007bl ). We start by mea- 
suring Ni, the number of galaxies within a cylinder of ra- 
dius 1 /i~ 1 Mpc and depth 120 h~ Mpc centred on every cen- 
tral galaxy in our catalogues (the depth mimics a redshift 
uncertainty of Az ~ 0.02). Then we discard those galax- 
ies whose cylinder overlaps with that of a central galaxy 
of a more massive cluster. This is equivalent to assuming 



that the luminosity of the central Brightest Cluster Galaxy 
(BCG) increases monotonically with mass, and then discard- 
ing as potential group centres those BCGs that are close 
to a brighter galaxy. After this cleaning procedure we use 
the galaxy counts around the remaining central galaxies to 
define a new "observed" cylinder radius R(N opt ), equal to 
the mean virial radius of clusters of the same Ni richness. 
Then, we repeat the counting and cleaning processes until 
we reach convergence. We keep all clusters down to a count 
of one galaxy (i.e. just the central BCG). We refer to this 
count as the optical richness N op t of the cluster. 

A serious systematic in optical cluster catalogues may 
be caused by misidentification of the BCG, and hence of the 
cente r of the corresponding dark matter halo (|Rozo et al.l 
l201lf ). This effect is referred to as 'miscentering'. In 
fact, 20% to 40% of the MaxBCG groups (depending on 
the cluster richness) suffer from this effect according to 
Ijohnston et alj (|2007T ). In order to mimic this in our anal- 
ysis, we have carried out our cluster identification pro- 
cedure after randomly displacing 30% of our candidate 
cluster centres according to a 2D Gaussian with a mean 
shi ft of 0.4fe~ 1 Mpc. Th is distribution of offsets is based 
on Ijohnston et al.l (|2007Tl , who applied the MaxBCG algo- 
rithm to mock cata logues built from the Hubble simulation 
(|Evrard et al.l |2002| ) . The functional form and parameters 
we use to describe the effect are also consistent with the dis- 
tribution of projected distances between the position of the 
dominant subhalo and that of the second most massive sub- 
halo within haloes of M XXL simulation (see also Fig. 2 of 
Iffilbert fc White! l201Ch . We caution that the miscentering 
fraction and the displacement parameters are uncertain and 
are sensitive to details of the clust er-finding algorithm. In 
iPlanck Collaboration et alj (1201 lal~) t he large X-ray cluster 
compilation of Piffaretti et al.l i 201 lh was matched to the 
maxBCG catalogue, finding a median offset between X-ray 
centre and BCG position of about 100 kpc for the 189 clus- 
ters in common; ~ 15% are offset by more than 400 kpc (J-B 
Melin, priva te communication) . Th is agrees reasonably with 
the result of Ijohnston et all (|2007l ). We will see below that 
few of our results are sensitive to miscentering because it 
generally causes cluster observables to be perturbed parallel 
to the scaling relations which link them. Since we retain a 
flag which notes which of our clusters have had their centres 
displaced, we can construct a cluster sample based on "true" 
centres simply by ignoring these objects. 

Our catalogue contains 594399 objects with JV op t above 
10, corresponding roughly to haloes with M200 > 4 x 
10 13 /i _1 M Q . There are 1988 objects with more than 100 
members, corresponding to M200 > 7 x lO 14 /i -1 M0. For 
each cluster we compute associated surrogate observables, 
the X-ray luminosity Lx, the weak-lensing mass Mi ons , and 
the tSZ flux Ysz, by integrating the corresponding 2D maps 
around the apparent (i.e. after 'miscentering') centre out to 
R(N pt). For each signal we subtract the contribution of the 
background, which we estimate using an annulus of radius 
1.5 x ~R(N opt ) < r < 2 x R(N opt ). Naturally, this is not the 
approach that one would follow for individual well-observed 
clusters, where one can directly identify the peak of the X- 
ray, tSZ or weak lensing signals and estimate an individual 
virial radius from a profile built around this centre. It is, 
however, closely analogous to the procedure followed when 
estimating scaling relations (mean values of Mi ens , Lx or Ysz 
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Figure 5. Virial masses (M200) for a simulated volume-limited sample of optically detected galaxy clusters at z = 0.25 as a function of 
our surrogate observables. Darker regions correspond to a larger number density of clusters. Explicitly, N op t corresponds to a maxBCG- 
like optical richness, Lx to projected X-ray luminosity, Ygz to projected tSZ signal and M\ cns to weak lcnsing mass. All these signals 
are integrated values within the virial radius, calculated as described in the text, using cluster centres that are randomly displaced 30% 
of the time from the true potential minima. Dot-dashed lines indicate the expected 'self-similar' scalings, while blue solid lines indicate 
the regressions (in log-log space) of halo mass against each surrogate. These have the slopes listed in each panel which differ from the 
self-similar expectation. Circles show the mean virial mass for a series of narrow logarithmic bins of each surrogate, while dashed red 
lines indicate the region containing the central 68% of the clusters in each bin. The fractional scatter in halo mass at given surrogate 
value is given in each panel. 



as a function of optical richness) by stacking large samples 
of optically selected clusters. 

In our analysis below, we consider two types of clus- 
ter sample that mimic catalogues from large-scale optical 
and X-ray surveys. By comparing results from these sam- 
ples we hope to assess the impact of the observational selec- 
tion method on deriv ed scaling relations (c.f. section 3.2). 
iKoester et al.l |2007al ) show that, to a good approximation 
and over a wide range of richness, their maxBCG catalogue 
can be considered to be selected over a fixed volume of about 
0.5 Gpc 3 . We thus take "optically selected" samples to be se- 
lected uniformly from the full MXXL volume, independently 
of their properties. 

Current large "representative" surveys of X-ray clusters 
are based primar ily on the Rosat All Sky Survey (RASS, 
IVoges et al.lll999l ) and on serendipitous discovery in fields 
observed for other reasons by the Rosat, Chandra and 



XMM-Newton satellites. The largest compilation s, such as 
the 1800 cluster MCXC of IPiffaretti et all (|201ll ). are built 
by combining subcatalogues each of which is effectively X- 
ray flux limited, and so surveys a substantially larger volume 
for X-ray bright clusters than for X-ray faint ones. For each 
subcatalogue, and so for the compilation as a whole, the 
volume surveyed scales approximately as L x because the 
apparent luminosity decreases as distance squared whereas 
the volume increases as distance cubed. This overrepresen- 
tation of luminous objects is known as Malmquist bias, and 
we incorporate it in our "X-ray flux limited" samples by re- 
taining all objects in the MXXL volume but weighting each 
by the 3/2 power of Lx- 
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Figure 6. Correlations among deviations of observables in galaxy clusters with mass in the range 4 X 1O 14 M < M200 < 1 X 10 15 M o . 
Data correspond to logarithmic deviations, i.e. Alog(s) = log(s) — (log(s)), where the mean is computed for clusters in narrow mass bins 
(A log M200 = 0.2). The intensity of the background 2D histogram is proportional to the number of haloes in the corresponding region 
of the plot, with a darker grey-scale indicating a larger number density of objects. Red circles correspond to the average y-value in bins 
along the x-axis. The linear correlation coefficient r for each pair of observables is given in the legend of each panel. 



3.2 Scalings with mass 

In Fig. [5]we present the relations between cluster virial mass 
M200 and each of our four surrogate observables for opti- 
cally detected clusters. Specifically, N opt corresponds to op- 
tical richness, Ysz and Lx to the projected tSZ and X-ray 
fluxes, and Mi ons to weak lensing mass, all integrated within 
the projected virial radius corresponding (in the mean) to 
its estimated richness and surrounding its apparent centre 
(i.e. after "miscentering" perturbations). The symbols in 
each panel correspond to the mean of the AL200 distribu- 
tion at given value of the surrogate observable and the red 
dashed lines contain the central 68% of this distribution. 



The straight blue lines show the linear fit to the individual 
cluster data (in logarithmic space) which minimises the rms 
residuals in the vertical direction. 

For comparison, green dot-dashed lines indicate the 
scaling expected naively given our assumptions about the re- 
lations between baryonic and dark matter properties. In the 
optical case, this comes directly from the HOD model used 
to build the galaxy catalogues (N opt oc Af 0,9 ), in the X-ray 
and tSZ cases from standard self-similar scaling (Lx oc M 4 ^ 3 
and Ysz oc M 5 ' 3 , respectively), and in the case of lensing it 
is direct proportionality (Mi cns oc M). In all four panels, 
the slope of the measured regression, shown in the legend, 
is similar but not identical to the expectation. These devia- 
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Figure 7. Correlations among deviations in Lx> ^SZ an( i ^lens 
for objects with optical richness in the range 44 < N opt < 50. 
The intensity of the background 2D histogram is proportional to 
the number of haloes in the corresponding region of the plot, with 
a darker grey-scale indicating a larger number density of objects. 
The linear correlation coefficient r for each pair of observables is 
given in the legend of each panel. 



tions can be explained by the relative impact for the different 
surrogates of internal halo structure and of contamination 
along the line-of-sight, as well as of miscentering. All these 
can depend systematically on halo mass. The largest dis- 
crepancy is found for the lensing surrogate, followed by the 
optical richness. The smallest are found for the X-ray and 
tSZ signals. This is consistent with the fact that the stronger 
the dependence of a surrogate on mass, the less sensitive it is 
to contamination and to other projection effects, since these 
are typically produced by less massive systems. 

The scatter in halo mass at given value of an observ- 
able can be roughly described by a log-normal distribution 
and depends weakly on the actual value of the observable 
for those we study here. It is indicated in each panel and 
ranges from 20 to 40%. The tSZ signal shows the least scat- 
ter and the lensing the largest. The values we find are con- 
sistent with previous studies, but note that our estimator s 
are not optimal fc.f. iMelin et al.ll2006l ; iRvkoff et alj|2012h . 
Note also that our catalogues do not include all possible 
sources of scatter, so even larger values may apply to real 
data. On the other hand, we consider that our results should 
yield a reliable upper limit on the size of uncorrelated pro- 
jection effects, given the very large box size of the MXXL. 
In particular, for tSZ our scatter estimates agree with those 
reported from full hydrodynamical simulation s of smaller 
volumes (|Kav et al.ll201ll ; Battaglia et al. 201lT); and for op- 
tical richness, w ith those directly inferred from the data for 
optical clusters (|Rozo et al.ll2009l ). This is a reassuring con- 
firmation that our assumptions are reasonable. 

An interesting corollary of the considerable dispersion 
in these relations is that it is unlikely that the cluster with 
the largest value for any particular observable will actually 
correspond to the most massive halo in the survey which 
identified it. We have explicitly checked that this perhaps 
counterintuitive situation does indeed hold. In the insets of 
Fig. [TJ we show the most extreme cluster in our simula- 
tion as identified by each of our four surrogate observables, 
i.e. the cluster with the largest tSZ signal, X-ray flux, grav- 
itational lensing signal and optical richness count. These 
clusters turn out all to be different. The cluster with the 
largest richness is, in fact, the one with the largest M200 at 
z = 0.25 and is notable also for the fact that it does not even 
have a well defined centre. This makes clear that consider- 
able care is needed to draw cosmological inferences from the 
observed properties of the most extreme cluster in any par- 
ticular survey. Statistically meaningful constraints can be 
obtained only with a complete and accurate treatment of 
the scatter in the mass-observable relation, including any 
possible dependence on cluster mass. 

An important point for our subsequent analysis is that 
while some sources of scatter affect primarily, or even exclu- 
sively, one specific observable, most affect several simulta- 
neously. For instance, at fixed M200 the HOD is expected to 
correlate with halo formation time because older and more 
relaxed haloes tend to have more dominant BCGs. Since 
formation time correlates strongly with concentration but 
only weakly with virial temperature, X-ray luminosity is also 
expected to increase with halo age, whereas integrated SZ- 
strength (and also lensing mass) should be age-independent. 
It has long been clear observationally that at given richness 
more regular and relaxed clusters (hence "older" clusters) 
do indeed have more dominant cD galaxies, and it is now 
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clear that they are also more X-ray luminous. In contrast, 
variations in baryon fraction are expected to affect Lx and 
Y~sz similarly, but to have little effect on Mi en s and to corre- 
late in a model-dependent and uncertain way with JV op t. In 
addition, orientation is expected to have little effect on the 
measured X-ray luminosity of a cluster but to produce corre- 
lated variations in its measured SZ flux, richness and lcnsing 
mass. Finally, misidentification of the centre and misestima- 
tion of the virial radius of a given cluster will induce vari- 
ations in all its observables. Generically, such effects imply 
that deviations from the various mass-observable relations 
are not independent. Rather, there is a non-zero covariance 
which reflects common sensitivities to halo structure, orien- 
tation, environment and foreground/background superposi- 
tion - surrogates which are similarly sensitive to these fac- 
tors are expected to exh ibit a high degree of correlation (see 
also IStanek et ai1l2010l ') . 

We quantify this effect in Fig. [6] which shows scatter 
plots of the deviations from the mean at given M200 in the 
logarithms of the values of observables for individual clus- 
ters. Here we include all clusters with 1 x 10 M@ > M200 > 
4 x 1O 14 M0. In each panel we give explicitly the Pearson 
correlation coefficient, r, which characterises the correlation 
between the deviations. 

The strongest correlation is that between the devia- 
tions in lensing mass and Ysz, presumably because they are 
similarly sensitive to cluster orientation, projection, miscen- 
tering and misestimation of ifeoo- The second strongest is 
between Lx and Ysz, the two quantities sensitive to our 
estimates of gas density and temperature. The weakest is 
between richness and Lx, perhaps because the X-ray lumi- 
nosity is dominated by the dominant central concentration 
of clusters while N opt is influenced substantially by orienta- 
tion and projection effects. The other three correlations are 
all of similar strength. 

While Fig. [6] illustrates the correlated scatter in observ- 
ables among clusters of given 'true' mass, the more relevant 
correlations for the effects discussed in the next section are 
those at fixed observed richness, N op t. These are shown as 
scatter plots for 44 < iV op t < 50 in Fig. [7] We have checked 
and found quite similar results for other richness ranges. 
The scatter in each observable is considerably larger at fixed 
7V op t than at fixed M200 and the correlations are substan- 
tially stronger in Fig. than in Fig. [6] reaching r = 0.85 for 
the particularly relevant case of Ysz versus Lx. 

Despite the different degrees of covariance among our 
surrogate observables, we note that we measure a positive 
correlation in all cases. This means that a cluster with an 
abnormally high signal in one surrogate is likely to have 
a high signal in the other three as well, especially those 
where the relevant correlation is strong. Hence, different ob- 
servables do not provide independent measurements of the 
true mass of any given cluster, and any analysis which as- 
sumes that they do is likely to be at least partially in error. 
Another implication is that a group of clusters selected us- 
ing of one of these observables will not form an unbiased 
sample of the underlying cluster population with respect to 
any of the other observables. As a result, both the mean 
scaling relations for such a sample and the scatter around 
these relations may differ from those for the full underlying 
cluster population. We will see, for example, that the mass- 
observable and observable-observable relations derived from 
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Figure 8. Thermal SZ signal Ygz as a function of optical richness 
N op t for a volume-limited sample of ~ 500, 000 clusters from the 
MXXL. The upper panel shows the mean relation and the lower 
panel the ratio of the scatter to the mean. In each panel the red 
curve shows the "intrinsic relation" where both the tSZ signal 
and richness are computed in 3D within a sphere of radius R200 
centred on the potential minimum. For the "observed" lines, both 
Ysz and 7V op t are computed in projection about these true cen- 
tres using the true -R200 an d using only uncontaminated clusters, 
defined to be those where at least 90% of the galaxies are mem- 
bers of a single FoF halo. The "+contaminated" lines refer to 
samples where this latter condition is eliminated and where the 
virial radius for each object is estimated iteratively from JV op t as 
is done for real data. Finally, the "+miscentred" curves refer to 
samples where the projected position of the centre has been offset 
for 30% of the clusters as described in the text. 

an optically selected cluster catalogue will differ in general 
from those derived from an X-ray selected catalogue. 

3.3 Systematic effects in measured scaling 
relations 

We now consider how the thermal SZ signal is related to 
optical richness for cluster samples selected in various ways. 
In Fig. [8] we display results where each MXXL cluster is as- 
signed equal weight in order to mimic volume-limited sam- 
ples like the optically selected maxBCG survey. For compar- 
ison, in Fig.|9]we show results where each cluster is weighted 
3/2 

by L x , mimicking cluster samples like the MCXC which is 
constructed from a number of X-ray surveys, each of which 
is effectively X-ray flux limited. The upper panels in these 
plots show mean Ygz for clusters of given N opt , while the 
lower panels show the rms scatter about these relations ex- 
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Figure 9. SZ signal Ygz as a function of optical richness N op t for 
a X-ray flux limited sample of ~ 500, 000 clusters from the MXXL. 
This is directly analogous to Fig. \E\ and was constructed in the 
same way except each cluster is weighted by the 3/2 power of its 
X-ray luminosity. The "observed" , "+contaminated" and "+mis- 
centered" lines refer to the same cluster sets as in Fig. [8] differing 
only because of this weighting. The red lines labelled "volume lim- 
ited" repeat the "+miscentered" results from Fig. [8] Malmquist 
effects substantially enhance the amplitude of the mean relation, 
reduce the scatter, and make the relation insensitive to superpo- 
sition and miscentering effects. 

pressed as a fraction of the mean signal. Lines in the upper 
panels thus represent the mean tSZ signals expected if clus- 
ters from optical or X-ray surveys are stacked as a function 
of their optical richness. 

The solid red lines in Fig. [H] show the "intrinsic" rela- 
tion which is obtained if both Ysz and N opt are calculated 
in 3D by integrating over a sphere centred on the poten- 
tial minimum and with radius ifeoo- The properties of real 
clusters are, of course, measured from 2D maps. The dotted 
lines labelled "observed" in both figures show the relations 
obtained for uncontaminated clusters when both Ysz and 
7V pt are integrated over a disk around the potential mini- 
mum with radius given by the true ifeoo ■ Here "uncontami- 
nated" means that at least 90% of the galaxies contributing 
to iVopt have to lie in a single FoF halo. Although, by defi- 
nition, our estimates of N opt and Ysz both increase for any 
individual object in going from 3D to 2D, they increase by 
similar amounts, with the result that the apparent relation 
does not change significantly over the full range of richness 
probed here. The scatter, on the other hand, is greatly in- 
creased for poor clusters, where it can double the intrinsic 
value, but is little affected for rich clusters. 



For the dashed blue lines labelled "-{-contaminated" in 
Figs. [8] and [9] we relax the requirement that the clusters 
be uncontaminated. We also use an observational estimate 
of -R200 for each cluster when estimating its richness and 
its tSZ signal, as described earlier in this section. Line-of- 
sight superposition effects then contribute to the values of 
all the observables. Because the superposed objects are usu- 
ally of lower mass (and hence lower temperature) than the 
main cluster, they typically inflate iVopt (which scales ap- 
proximately as M ' 9 ) more than they do Ysz (which scales 
approximately as M 5 ^ 3 ). As a result, the "-(-contaminated" 
relations lie slightly to the right of the "uncontaminated" 
ones. The effect is smaller in Fig. [9] than in Fig. [8] The 
scatter is further increased by contamination in the volume- 
limited case but is relatively little affected for flux-limited 
samples. 

The final observational systematic we study is miscen- 
tering. The dot-dashed purple curves labelled "-fmiscenter- 
ing" in Figs. [8]and[9]show the relations found when the cen- 
tres assigned to a random 30% of the clusters are offset from 
their potential minima as described above. This results in a 
surprisingly small shift in the mean relation in both cases. 
Again, this is because such offsets induce shifts in the esti- 
mated values of iV op t and Ysz that are largely parallel to the 
mean relation. 

These "+miscentering" curves give our most realistic 
estimate of the relations expected for real clusters in the 
two cases. We repeat these (dot-dashed purple) curves from 
Fig. [8] as solid red curves in Fig. [9] in order to emphasise the 
most important result of this section. The average Ysz signal 
for X-ray flux-limited samples is boosted by a factor of 3.5 in 
the low richness tail and by a factor of 1.25 at high richness, 
relative to the volume-limited case. This is a consequence 
of the strong correlation between _Lx and Ysz at fixed N opt 
which is visible in the top panel of Fig. In samples selected 
above a limiting X-ray flux, clusters of given N op t which are 
X-ray underluminous are down-weighted and these tend also 
to be the objects with the smallest Ysz signals. This is a 
manifestation of Malmquist bias. 

A second important consequence of X-ray selection is 
that the scatter about the mean Ysz-A^pt relation is greatly 
reduced. The lower panel of Fig. [5] shows that we predict 
it to be only about half that for a complete volume-limited 
sample of clusters. Although much of this difference is due to 
a reduced sensitivity to contamination and miscentering, the 
scatter is lower even than that shown in Fig.[8]for uncontam- 
inated and properly centred clusters. If not corrected, this 
could lead to over-optimistic estimates of the performance 
of mass estimators when tested on X-ray selected cluster 
samples. Such samples are clearly more homogeneous in in- 
ternal structure at given mass than volume- limited samples. 
We expect biases of this kind to be present in any cluster 
catalogue selected according to a specific observable, and 
they must be corrected in order to infer correctly, for exam- 
ple, the volume abundance of clusters as a function of M200, 
the quantity normally used to draw cosmological conclusions 
from the cluster population. Robust constraints require not 
only that the mean transformation from observable to mass 
be determined accurately and without bias, but also that the 
scatter between these quantities be known to high precision. 

Fig. [10] is analogous to the previous two figures but now 
focuses on the relation between gravitational lensing mass 
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effect in X-ray flux limited cluster samples because at given 
richness more luminous clusters are usually more massive. 

The various biases examined here depend both on the 
intrinsic properties of the cluster population and on the 
specific technique used to measure each observable. For in- 
stance, the smaller the scatter between X-ray luminosity and 
optical richness, the smaller the impact of Malmquist bias 
on the amplitude of the Ysz-N op t relation. In the limiting 
case of no scatter the relations for flux-limited and volume- 
limited samples would be identical. The same would be true 
if there were no correlation between Lx and Ysz at fixed 
iVopt- The substantial effects found above are due to the 
strong correlations we predict. The details of the observa- 
tional procedures matter because they can enhance or sup- 
press the impact of different aspects of the data. For in- 
stance, an optimal filter based on expected cluster profiles 
will be more sensitive to miscentering than a top-hat filter 
of larg e size; the shear profile scheme which Uohnston et al.l 
120071 ) used to estimate lensing masses may be yet differently 
sensitive. Such effects must be taken into account properly if 
any particular survey is to place reliable constraints on cos- 
mological parameters. This also applies to inferences from 
the properties of extreme objects. Cosmological inferences 
require a full understanding of the scatter in the observable- 
mass relation and an accurate knowledge of the selection 
function, otherwise one may arrive at seriously erroneous 
conclusions. 



Figure 10. Mean lensing mass as a function of optical richness 
for clusters in the MXXL simulation. "Lensing mass" is defined as 
the mass within a properly centred sphere of radius i?200 (-^200 > 
the red solid curve labelled "intrinsic" ) or projected within a cir- 
cle of radius given by an iVopt-based estimate of R200 and centred 
cither on the true potential minimum ( "observed" and "-(-contam- 
inated" cases) or on an offset centre in 30% of the cases ( "-(-mis- 
centred"). N op t is defined within these same regions. The upper 
panel shows the mean relations while the lower one gives the 
fractional scatter about these relations. The "X-ray flux limited" 

curves refer to the "+miscentred" case with haloes weighted by 
,3/2 



and richness. As was the case for Ysz-Aopt, we find that 
the slopes of the mean relations are similar in all the vol- 
ume limited cases, and that the effects of co ntamination and 
miscen tering are quite modest. In contrast to ljohnston et aD 
(2007). who find that miscentering decreases the normalisa- 
tion of the Miens-^opt relation by 15-40% (see their Tables 
3 and 8), it barely changes the mean relation in our data. 
This is because of the strong correlation between our esti- 
mators of optical richness and lensing mass. If the center of 
a halo is misidentified, both JV op t and Mi ona are underesti- 
mated, and the change is, on average, almos t parallel to the 
mean relation. The stronger effect seen by Ijohnston et alj 
{2OO7]) may reflect their different estimator for the lensing 
mass, or possibly a failure to account consistently for the 
implied change in richness. We note that this richness re- 
duction implies that at any given Ar op t, fewer than 30% of 
clusters will actually be miscentered, both because cluster 
abundances increase steeply with decreasing richness, and 
because poor clusters are more likely to be rejected by our 
algorithm in favour of richer overlapping systems. Note also 
that the mean lensing mass is biased high by the Malmquist 



4 SCALING RELATIONS FOR PLANCK 
CLUSTERS 

In the last section, we showed that cluster scaling rela- 
tions are strongly and systematically affected by the way 
in which cluster samples are selected. This is a consequence 
of the substantial scatter in mass-observable and observable- 
observable relations, and the fact that deviations of differ- 
ent observables from the mean relations are strongly corre- 
lated because of common sensitivities to cluster structure, 
orientation, environment, and line-of-sight projection. The 
resulting distortion of scaling relations depends on how clus- 
ters are detected and their observables measured, so precise 
correction requires detailed modelling of each individual ex- 
periment. Our catalogues based on surrogate observables are 
nevertheless realistic enough to examine whether effects of 
this kind might explain some apparent inconsistencies re- 
cently highlighted by the Planck Collaboration. 

There are three pieces to the puzzle presented by the 
Planck collaboration. First, the mean tSZ signal for stacks 
of maxBCG clusters of given optical richness is about half 
of that predicted by scaling relations derived from the much 
smaller REXCESS sample for which in dividual X-ray pro- 
files are available jBohringer et al.ll2007l ) . The hot gas struc- 
ture of the REXCESS sample is well described by the al- 
most self-similar scaling of a "universal" pressure profile 
which resembles that predi cted by hydrodyna mical simu- 
lations of cluster formation |Arnaud et al.ll2010l ). The REX- 
CESS scaling relations do, however, agree well with the mean 
Ysz measured when cluste rs from the large MCXC compila- 
tionj Piffarctti ct al. 201 1|) are stac ked as a function of Lx 
|Planck Collaboration et al.l [20lTbh . Note that comparison 
with the stacked Ysz-N opt relation for maxBCG sample re- 
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Figure 11. Comparison of scaling relations between observed and simulated clusters. The left-hand panel shows X-ray luminosity as 
a function of optica l richness. The red circles with error bars show the mean luminosity of stacks of MaxBCG clusters as measured by 
IRvkoff et al.l d2008ah . while blue points indicate individual maxBCG clusters which are also in the MCXC compilation of lPiffaretti et al] 
j201lh . The right-hand panel shows the tSZ signal as a function of optical richness. Here, red circles indicate mean values for stacks of the 
entire maxBCG catalogue, while blue circles are for stacks of maxBCG clusters that arc also in the MCXC ( Planck Collaboration et al. 
l201ld ). All observational data correspond to the equivalent of measurements within -R500 and at z = 0. Predictions from the MXXL 
simulation are shown by red lines for a volume-limited sample and by blue lines for a X-ray flux-limited sample. Note that the normali- 
sations of our Lx an d Ygz proxies are set by shifting the red lines vertically to give the best possible fit to the red circles in these two 
plots. 



quires an estimate of cluster mass as a fu nction of N opt , 



for 



which lPlanck Collaboration et all J201ld) used weak lensin 



results for stacke d maxBCG clusters ( Johnston et al.ll2007l ; 
iRozo et alj|2009l ). 

The second piece of the puzzle is that if the maxBCG 
sample is restricted to clusters which are also in the MCXC, 
the stacked Ysz-N opt relation lies well above (a factor of 2 
at iVopt = 50) that for the sample as a whole and appears 
consistent with the REXCESS scaling relations. The third 
and final piece is that the relation between mean Ysz and 
mean Lx signals for stacks of the full maxBCG binned by 
A^opt is consistent both with that between mean (stacked) 
Ysz and Lx in the MCXC and with that between Ysz and 
Lx for individual clusters in the REXCESS sample. All these 
facts led the Planck Collaboration to speculate that a sub- 
set of optically detected clusters might have very weak tSZ 
and X-ray signals, presumably because they contain many 
galaxies but rather little hot gas. 

We now examine these issues using suitably selected 
cluster catalogues from the MXXL. To mimic the full 
maxBCG catalogue, we use volume-limited samples (c.f. 
iKoester et al.l l2007aT ) and we measure a "maxBCG-like" 
richness for each as detailed in Section 3.1. To mimic the 
MCXC, we will use "X-ray fl ux limited" samples which 
weight each object by L x , since iPiffaretti et~aH (|201ll ) built 
the MCXC by combining a number of X-ray surveys, most 
of which are effectively flux-limited. We are also interested 
in the overlap between these two observational samples. At 
high X-ray luminosity clusters are detected in the RASS to 
distances beyond the limit of the maxBCG catalogue, so 
that the combined sample is effectively volume-limited. At 



low X-ray luminosity, on the other hand, the maxBCG limit 
is well beyond the distance at which RASS can detect clus- 
ters and the overlap sample is flux-limited. As we will see 
below, the effects of this change in sample selection are di- 
rectly visible in the behaviour of the stacked Ysz signal of 
the overlap sample. 

We begin by considering the mean X-ray luminosity and 
mean tSZ signal as functions of optical richness for MaxBCG 
clusters. These data are displayed in Fig. [11] as red cir- 
cles. In the l eft pan el, we show X-ray measurements from 
IRvkoff et al.l |2008al ) which correspond to the average X-ray 
luminosity within -R500 for stacks of RASS maps centered 
on MaxBCG clusters. We convert their luminosities, origi- 
nally measured within -R200, to ones measured within -R500 
by multiplying them by 0.91. We also scale these values to 
their z = equivalent by assuming a self similar scaling 
of the luminosities. In the right panel we show the average 
tSZ signal in Planck maps within -R500 for the same richness- 
binned M axBCG clus t ers. N ote that the mass-richness rela- 
tion from IRozo et al.l (|2009h was used to set the size of the 
matched filter when measuring the Ysz signals, but that this 
has little impact on the result and is unimportant for our 
discussion here. 

We now compare these results to our simulation. Red 
lines show the mean relation predicted by "observational" 
catalogues which include both contamination and miscen- 
tering. They have been shifted vertically to fit the maxBCG 
data as well as possible, thus determining the scaling coeffi- 
cients used consistently in our Lx and Ysz proxies through- 
out this paper. After this normalisation, the Ysz — N200 and 
Lx — N200 relations are predicted remarkably well. The lat- 
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ter relation in the MXXL is slightly shallower than observed. 
This reflects the well known fact that the observed Lx-Tx 
relation is steeper than that expected from the self-similar 
model on which we built our surrogates, apparently indicat- 
ing a systematic variation in gas fraction or concentration 
with cluster mass. The agreement nevertheless appears good 
enough to validate our approximate modelling of gas physics, 
suggesting that it is adequate to capture the statistical cor- 
relations underlying the influence of selection bias on cluster 
scaling relations. 

The second set of observations in Fig. [11] refer to the 
set of 189 maxBCG clusters which are also in the MCXC. 
Blue dots in the left-hand panel show Lx as a function of 
-Nopt for the individual clusters, whereas blue squares in the 
both panels give average values for stacks of these clusters 
around their maxBCG centres, using the same N OP t bins 
as for the full maxBCG sample jPlanck Collaboration et alj 
l201lj >. In both cases we plot the z = equivalent value, 
using the measured redshift for each cluster and assuming 
a self-similar scaling of the signals. This allows a consistent 
comparison of all datasets. 

As noted by the Planck Collaboration, stacked tSZ 
fluxes are systematically larger for MCXC clusters than for 
the full MaxBCG sample, except possibly for the richest 
systems. The left panel ofFig. [TT] indicates that their X- 
ray luminosities are also systematically higher, again with 
the possible exception of the richest clusters. Blue dashed 
lines in both panels indicate the mean relations we predict 
for X-ray flux limited samples. The Malmquist bias offset 
from the volume-limited relation explains part of the differ- 
ence between the blue squares and the red circles in the left 
panel, and it explains almost completely the discrepancy in 
the right panel. As discussed in the last section, the latter 
is caused by a strong correlation between the deviations of 
individual clusters from the mean Lx-N op t and Ysz-A op t re- 
lations. This causes Malmquist bias to propagate from X-ray 
selection into the Ysz- N op t relation^. This suggests that the 
correlated scatter between Lx and Ysz at given N op t is well 
represented in our model, but also that there are sources 
of scatter which only affect Lx that are not accounted for 
in our analysis which would increase the difference between 
X-ray scaling relations derived from volume and X-ray flux 
limited samples without affecting the Ysz signal. 

Finally, we reiterate that both relations have consider- 
able scatter. At given N op t our model indicates that the frac- 
tional uncertainty in Lx and Ysz for volume-limited samples 
is about 40% for N opt ~ 200 and rises to 130% for iV opt ~ 10. 
These numbers are broadly consistent with the intrin sic scat- 
ter reported by IPlanck Collaboration et al.l (|201ld ) for the 
SZ measurement of MaxBCG clusters. The corresponding 
fractional scatters are about 40% for both Lx and Ysz in 
our X-ray flux limited samples. 

We now move to another observable scaling relation re- 
ported by PLANCK: that between Ysz and Lx- We display 
simulation results for both volume-limited and X-ray flux 
limited catalogues in Fig. 1121 Clusters were stacked as a 
function of 7V pt and then the mean Ysz of each stack was 
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J iRvkoff et al. | i2008rj) give extensive discussion of various bias 
effects when the maxBCG catalogue is combined with X-ray clus- 
ter surveys. 



Figure 12. Average tSZ flux as a function of average X-ray lumi- 
nosity for clusters stacked according to optical richness. Red solid 
and blue dashed lines give MXXL results for our volume-limited 
and X-ray flux limited samples, respectively and refer to the same 
set of Aopt bins. Both Lx and Ygz are larger for the flux-limited 
sample, but the shift is almost parallel to the mean relation so 
that Malmquist bias has little vis i ble effe ct. Red circles are taken 
from IPlanck Col laboratio n et al. | bond) and refer to maxBCG 
clusters, whilst blue squares are the results for the subset of clus- 
ters present in the MCXC catalogue. In both cases the data are 
binned according to N op t- Error bars indicate bootstrap uncer- 
tainties in the mean Lx and Ysz signals in each richness bin. The 
gr een dotted line shows the prediction s of the X-ray model build 
by IPlanck Collaboration et al] l|201ld l. 



plotted against its mean Lx. Although both Ysz and Lx 
are substantially larger at given jV opt in flux-limited stacks, 
the two relations are almost identical. Malmquist bias has 
a negligible impact because the shifts are almost parallel to 
the mean relation. 

The same relation can be constructed for the (almost) 
volume-limited maxBCG catalogue by stacking RASS and 
Planck data fo r clusters binned by N opt . This ex ercise was 
carried out by IPlanck Collaboration et al.l (|201ld l and their 
result is overplotted in Fig. 1121 We also include, as a green 
dot ted line, the relation predict e d by t he X-ray model built 
by IPlanck Collaboration et al.l (|201ld ) based on observa- 
tions of REXCESS clusters and an observational N op t-Mi ens 
calibration. Despite the relations being built from different 
samples, they are very similar to each other and to our 
MXXL predictions, showing not only that the Ysz-Lx re- 
lation is insensitive to details of sample definition, but also 
that our proxies continue to represent the observables well 
in this plane. 

An Ysz-Lx relation for the effectively flux-limited 
MCXC-maxBCG overlap is also shown in Fig. [12] For the 
richest optical bins, it is compatible with all other relations, 
but the Lx signal is considerably larger for poor clusters. 
The relatively small number of systems per bin leaves room 
for this to be a statistical fluctuation (and indeed, we find 
that the median per bin agrees much better with the ex- 
pected relation). However, if real, it indicates that the dif- 
ferences in Ysz and Lx at fixed iV op t cannot be explained 
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by the exist ence of a subpopulation devoid of hot g Els 
suggested bv lPlanck Collaboration et al.l (|201ld ) ■ since this 
would move the points in Fig. [12] along a line of slope unity 
and so could not bring them back onto the relation for the 
full maxBCG sample. On the other hand, all data can be 
simultaneously explained by Malmquist bias, together with 
an extra source of scatter affecting only Lx (e.g. the presence 
or absence of cool cores in the gas distributions). Note that 
this requires correlated scatter in observables of the kind il- 
lustrated in Fig. [6] and Fig. [7| and does not work in simpler 
treatm ents of the same puzzle, as that of lBiesiadzinski et alj 
|2012l ) which does not fully incorporate all such correlations. 

Given that our simulation appears to reproduce consis- 
tently all the measured observable- observable relations for 
the maxBCG and MCXC catalogues, it is interesting to dis- 
cuss why the Planck collaboration's modelling gave a pre- 
diction which was inconsistent with the measured Ysz-N opt 
relation for the maxBCG sample. The answer appears to 
lie in the way cluster masses were used to relate the popu- 
lation identified in X-ray surveys to the maxBCG clusters. 
Since values of N op t are not available for most of the clus- 
ters in bright and well studied X-ray samples, this richness 
measure must be calibrated against another observable if 
the hot gas properties of such samples are to be used to 
predict the mean Ysz signal from m axBCG stack s . X-ray 
luminosity could have been used, since iRvkoff et ail |2008br ) 
provide mean Lx values as a function of N opt , and the 
Ysz-Lx relation is not only theoretically robust but also 
observationally well determined (see Fig. I12[) . This route 
(N opt — ?> Lx — > Ysz) predicts mean Ysz values for maxBCG 
stacks which agree well with those measured directly. How- 
ever, the Planck collaboration decided instead to follow a 
different route (N op t — ► M500 — > Ysz) using the mean weak 
lensing masses meas ured as a functi o n of iV pt for stacked 
maxB CG clusters by I Johnston et all (|2007f ) and lRozo et alj 
|2009l ). These studies made substantial upward corrections 
to the directly measured masses to account for the effects of 
line-of-sight contamination and miscentering, but they failed 
to make consistent corrections to ./Vopt, thus ignoring the cor- 
related deviations of individual clusters from the mean N op t- 
M200 and Mi cns -M2oo relations visible in the lower left panel 
of Fig[Sl These correlations ensure that the apparent N op t- 
Afiens relation is almost independent of misce ntering. In fact, 
if the original "raw" masses obtained by Ijohnston et alj 
(|2007l ) are used instead of the "corrected" masses to pre- 
dict mean Ysz as a function of N op t, the discrepancy with 
the directly measured values almost disappears. 



5 CONCLUSIONS 

Throughout this paper, and in particular in the previous sec- 
tion, we have emphasised the importance of understanding 
how survey methods influence the scaling relations measured 
in the galaxy cluster catalogues they produce. We have ar- 
gued that this is crucial both for proper statistical analysis 
of the physical properties of the cluster population and for 
deriving meaningful cosmological constraints from the esti- 
mated masses of the extreme clusters identified in any given 
survey. 

In order to illustrate these points, we have used the dark 
matter distribution in the MXXL simulation, the largest 



high-resolution cosmological calculation to date, to con- 
struct sky maps from which clusters can be catalogued us- 
ing proxies for four different observables: optical richness as 
measured in deep photometric redshift surveys, X-ray lumi- 
nosity, thermal Sunyaev-Zeldovich (tSZ) signal, and weak 
lensing mass. Although our treatment of these observables 
is necessarily simplified, it is sufficient to explore the scatter 
in the relation between each observable and cluster mass, 
as well as the correlations between the scatter for differ- 
ent observables caused by common sensitivities to the in- 
ternal structure, orientation, environment and background 
contamination of clusters. This is essential to understand the 
systematic biases imposed by specific observational strate- 
gies for detecting clusters and measuring their properties. 

We employed these catalogues to show that there are 
a number of effects that systematically alter the slope, am- 
plitude and scatter of scaling relations among the observ- 
ables. Structural complexities, orientation variations, super- 
position both of surrounding large-scale structure and of un- 
related foreground/background objects, and miscentering all 
increase the scatter in the Ysz, Lx, M\ ens and JV op t signals 
for given cluster mass. Relations between these observables, 
for example, the M\ ens -N pt or Ysz-Lx relations, can, how- 
ever, be much less affected because clusters scatter roughly 
parallel to the mean relation. In addition, Malmquist effects 
in flux-limited surveys not only bias the amplitude and re- 
duce the scatter in the mass-observable relation for the ob- 
servable used to select the sample, but also in those for other 
observables which have correlated scatter. The strength of 
such effects depends substantially on survey strategy and on 
the operational definition of the observables. 

Ignoring these bias effects can lead to serious difficul- 
ties in interpreting cluster data. As an example, we have 
considered a discrepancy recently highlighted by the Planck 
collaboration which concerns the mean tSZ and X-ray sig- 
nals measured for stacks of clusters identified from optical 
and X-ray surveys. Both signals are lower for optically se- 
lected clusters than predicted for their weak lensing esti- 
mated masses by a model which fits both individually ob- 
served and stacked X-ray-selected clusters. Our results sug- 
gest that the data are nevertheless in good agreement with 
predictions for a concordance ACDM universe, even if the 
gas properties of clusters are assumed to scale in a simple 
self-similar way with cluster mass. The discrepancy appears 
to reflect Malmquist bias propagating from the X-ray lumi- 
nosities to the tSZ signal through covariance in their scatter 
at fixed cluster mass. Malmquist bias has rather little im- 
pact on the mean Ysz — Lx relation, since clusters scatter 
almost along it. The discrepancy appears to have been ex- 
acerbated by applying a substantial miscentering correction 
to the mean Mi ons for the stacked clusters without applying 
a corresponding correction to the mean values of the other 
relevant observables. Our model suggests that together these 
effects may resolve the apparent puzzle. 

Although our analysis appears to explain the discrep- 
ancy both qualitatively and quantitatively, our explanation 
should still be regarded as provisional because of the detailed 
dependence of the effects on how the observables were ob- 
tained from the observational data. A firmer conclusion can 
only be reached through considerably more detailed mod- 
elling of the particular surveys involved. This should address 
not only survey design and cluster identification issues, but 
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also the specific algorithms (matched filters, etc.) used to 
measure the observables. Additional uncertainty comes from 
our schematic treatment of the baryonic physics, which un- 
doubtedly misses important aspects of the relation between 
the visible material and the underlying mass. It is never- 
theless clear that precision cosmology with clusters will be 
impossible without carefully tailored surveys with calibra- 
tion strategies that fully account for the multidimensional 
scatter between all the relevant observables and the fiducial 
cluster mass. Furthermore, linking the observations to the 
underlying cosmological model will require simulations that 
model all these statistical and astrophysical aspects to the 
required level of precision. Even with its limited treatment 
of the relevant astrophysics, the remarkable size and statis- 
tical power of the MXXL gives a foretaste of what should 
be possible in the future. 
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